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Abstract 

I perform lattice Monte Carlo studies of universal four-component fermion systems in one spatial 
dimension. Continuum few-body observables (i.e., ground state energies and integrated contact 
densities) are determined for both unpolarized and polarized systems of up to eight fermions 
confined to a harmonic trap. Estimates of the continuum energies for four and five trapped fermions 
show agreement with exact analytic calculations to within approximately one percent statistical 
uncertainties. Continuum many-body observables are determined for unpolarized systems of up 
to 88 fermions confined to a finite box, and 56 fermions confined to a harmonic trap. Results 
are reported for universal quantities such as the Bertsch parameter, defined as the energy of the 
untrapped many-body system in units of the corresponding free-gas energy, and its subleading 
correction at large but finite scattering length. Two independent estimates of these quantities 
are obtained from thermodynamic limit extrapolations of continuum extrapolated observables. A 
third estimate of the Bertsch parameter is obtained by combining estimates of the untrapped and 
trapped integrated contact densities with additional theoretical input from a calculation based 
on Thomas-Fermi theory. All estimates of the Bertsch parameter and its subleading correction 
are found to be consistent to within approximately one percent statistical uncertainties. Finally, 
the continuum restoration of virial theorems is confirmed for both few- and many-body systems 
confined to a trap. 
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I. INTRODUCTION 



Universal Fermi gases have garnered wide-spread attention following their recent realiza- 
tion in ultra-cold atom experiments. Perhaps the most physically interesting example of 
such a system is the unitary Fermi gas in three spatial dimensions. This system comprises 
a dilute mixture of spin-1/2 fermions interacting via a short-range interparticle potential 
tuned to produce an infinite two-particle scattering length. It is universal in the sense that 
the physical properties of the system are independent of the short-distance character of the 
interaction. As such, unitary fermions are not only relevant for ultra-cold atoms, but also 
for a variety of physical systems considered in other disciplines as well. 

Early on, the unitary Fermi gas was proposed as an idealized model for describing dilute 
neutron matter in neutron stars [Ij. The system had later been realized experimentally in 
ultra-cold cold atom experiments by exploiting properties of a Feschbach resonance [2H6]. 
More recently, unitary fermions have gained theoretical interest as an example of a nonrela- 
tivistic conformal field theory. A greater experimental and theoretical understanding of this 
strongly interacting and nonperturbative system has revealed that despite their simplicity, 
unitary fermions possess many rich and rather surprising physical properties. 

The unitary Fermi gas is simple in the sense that the system is characterized by a single 
physical scale, its density p. From purely dimensional considerations, the energy density of 
the interacting system must be given by 

£{p) = e^o(p) , (1) 

where So{p) is the energy density for free fermions, and the dimensionless constant of pro- 
portionality ^ is a nonperturbative universal quantity known as the Bertsch parameter [7j . A 
calculation based on density-functional theory [8] suggests that for the same system confined 
to a harmonic trap, the energy is given by 

i?-(Q) = v/ei?o-(Q) , (2) 

where Eq'^'^{Q) is the energy of the corresponding noninteracting system in the limit of 
large fermion number Q. The parameter ^ appearing as the square of the constant of 
proportionality in Eq. [2] is the same as that in Eq. [ij This nontrivial relation, including 
subleading corrections due to a finite fermion number, was independently confirmed from a 
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calculation based on a general coordinate invariant effective field theory description of the 
system piTU]. 

During the past decade, substantial effort has been devoted to determining the numerical 
value of the Bertsch parameter using a variety of analytical, numerical and experimental 
means (see, e.g., pLlJ and references therein). Although analytical calculations have become 
quite sophisticated, they generically possessed unquantifiable systematic errors. This is due 
in part to the lack of any small dimensionless parameters which might allow a perturbative 
expansion. Presently the most accurately determined value for the Bertsch parameter based 
on experiment is ^ = 0.376(4) [12]. The current best estimate from lattice Monte Carlo 
studies of an untrapped Fermi gas have yielded C, = 0.372(5) [T3] . 

More recently, interest has turned toward lower- and mixed-dimensional Fermi gases, 
which may in principle be realized experimentally using optical lattices. Although perhaps 
more difficult to realize experimentally, it has been demonstrated that such systems can 
exhibit universal properties of their own [14J. The focus of this study is a nonrelativistic 
four-component Fermi gas in one spatial dimension, with an attractive short-range four- 
body interaction exactly tuned to produce a zero-energy four-body bound state. In this 
limit, the one-dimensional system becomes universal in the sense that properties of the 
system are independent of the details of the interparticle interaction. Furthermore, the 
theory exhibits zero-temperature few- and many-body properties which are in many ways 
qualitatively identical to those of the three-dimensional spin-1/2 unitary Fermi gas 

As was the case for unitary fermions in three dimension, the four-component one- 
dimensional unitary Fermi gas is characterized by only a single scale, the density of the 
system. By dimensional analysis, the energy density of the untrapped system must obey 
Eq. [T], where the free gas energy density is given by 

Similarly, a calculation based on Thomas-Fermi theory along the lines of fSJ shows that the 
one-dimensional Fermi gas confined to a harmonic trap obeys Eq. |2] at unitarity, where 

E^iQ) = Iq'^ , (4) 
and 00 is the characteristic frequency of the harmonic potentiaj^ (see Appendix |A] for details). 
The Bertsch parameter ^ associated with the four-component Fermi gas need not be the same 

^ Throughout I work in units where h — 1. 
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as that of the two-component Fermi gas, and at present there are no theoretical arguments to 
support that they are. Numerical evidence reported in a companion paper, however, suggests 
that the Bertsch parameter for these two systems are in fact equal to within approximately 
one percent statistical uncertainties [T6] . 

Moving away from the unitary limit, the energy of the untrapped one-dimensional unitary 
Fermi gas, normalized by the the free-gas energy, can be described by a universal Bertsch 
function, which depends only on the dimensionless quantity y = {kpa)~^. Here, 

kp = np/A is the Fermi momentum for free fermions and a is the scattering length, a 
length scale which characterizes the long distance character of the four-fermion interaction 
potential. In the case of the four-component Fermi gas, the Bertsch function is shown to 
behave as [15] 



_6_ 

e - Cy \y\ « 1 (5) 



1 + ^ + • • • y « -1 



4i 

where the Bertsch parameter ^ was previously discussed, and the slope ( at unitarity is 
another unknown nonperturbative constant. The slope of the Bertsch function in the unitary 
limit can be expressed as 



where the quantity 



C = -{Airm) 



(7) 



da 

is known as the contact density, evaluated in the unitary limit. The contact density is a 
well-defined physical quantity for both few- and many-body systems, at zero- and finite- 
temperature, and even away from unitarity. It furthermore plays a fundamental role in 
various universal (Tan) relations [T7H2T] . 

The Hamiltonian for the one-dimensional universal Fermi gas is invariant under symmetry 
transformations generated by the Schrodinger algebra [221 123] • The theoretical implications 
of these symmetries have been explored in great detail [211 [25], and include: 

1. an operator-state correspondence, which relates the scaling dimensions of primary 
operators in free space to the energy levels (in units of the trap frequency) of the 
system confined to a harmonic trap; 



2. 



virial theorems for trapped systems, which relate the expectation value of the potential 
energy operator to the expectation value of the untrapped Hamiltonian (as well as their 
powers) ; 



3. a tower of states in the trapped system which are separated by 2uj (i.e., breathing 
modes) . 

These properties are generic, holding for both the few- and many-body systems, and may 
be explored numerically via Monte Carlo simulations. 

In this paper, I perform detailed numerical studies of few- and many-body systems con- 
fined to a harmonic trap and finite box, and present continuum limit extrapolated results for 
their energies and integrated contact densities. In Sec. |TT| I summarize the salient features 
of the lattice description used for this study, discuss the role of parameters appearing in 



the lattice theory, parameter tuning and the continuum limit. In Sec. |III[ I introduce a 
world-line representation for the partition function, and present explicit definitions for the 
physical observables of interest, including the energy of the system and integrated contact 



density (contact). In Sec. IV, I present details pertaining to the algorithm used to simulate 
the theory, the simulation parameters considered, and the generation of statistical ensem- 
bles. In Sec. |V], I discuss the analysis of data, including details relating to how continuum, 
infinite volume, and thermodynamic limit extrapolations were performed. In this section 
I present continuum extrapolated estimates for few- and many-body observables confined 
to a trap and a box, and for the case of many-body systems, present thermodynamic limit 
estimates for the energies and integrated contact densities. 

The many-body results for the Bertsch parameter presented in Sec. |V] were originally 
reported in [16j. In this paper, I go into greater detail regarding the analysis of those 
results, as well as present new results for the Bertsch parameter and integrated contact 
densities using ensembles of increased size. From the latter estimates, I determine the 
universal parameter as well as make a third determination of the Bertsch parameter using 



additional input from a density-functional theory calculation. In Sec. VI, I summarize the 
results of this study and provide some concluding remarks. Finally in Appendix |A| I derive 
the dependence of the trapped many-body energy on ^ and C, for a one-dimensional Fermi 
gas in the unitary regime, and provide a confirmation of Eq. |2} In Appendix [B} I provide 
plots of all the continuum limit extrapolation results for few- and many-body observables of 



interest. 

II. THEORY 

The starting point for this study is an effective field theory for nonrelativistic fermions 
interacting via an attractive four-body contact interaction. The continuum Lagrangian for 
the theory defined in two-dimensional Euclidean space-time is given by: 



where iIj„{t,x) is a four-component Grassmann- valued spinor with spin components labeled 
by the index cr = (a, b, c, d) and space-time coordinates labeled by the coordinate pair (r, x), 
m is the fermion mass, and (7 is a coupling associated with the four-body interaction. I 
consider the theory defined in the presence of a spin-independent external potential f (x), 
given by the two choices: 



where the parameter k denotes the oscillator spring constant associated with a trapping 
potential. 

I consider the theory at a finite temperature T = employing anti-periodic boundary 
conditions in the time direction, and consider a system with open boundary conditions in 
the space direction. The choice for the latter is arbitrary, with different choices leading to 
different finite volume artifacts which are ultimately removed in the infinite volume limit. 
The continuum theory is discretized on an x {2Ns + 1) rectangular lattice with lattice 
sites labeled by the integer coordinate pair n = (n,-, ng) for rir G [0, Nr) and G [— A^s, Ng]. 
In lieu of continuous fields, one considers fermion fields ipn (and their Hermitian conjugates 
ipl^ defined only at the sites of the lattice. The continuum operators appearing in Eq. [sjare 
then defined on the lattice using conventional finite difference discretizations, following |2S] '■ 







T 



(10) 
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where br (bg) is the temporal (spatial) lattice spacing with r = fe^-n^ {x = bgns), and 
(Cs) is a unit vector pointing in the time (space) direction. The lattice-discretized external 
potential is given by = ^{bsnsY- The physical spatial extent of the lattice is given by 
L = bs{2Ns + 1), and the temporal extent (i.e., inverse temperature) is given by /3 = brNr- 
At infinite volume and at zero temperature, the ground state energy E of the untrapped 
four-body system can be analytically related to the four-body coupling by exact diago- 
nalization of the four-body transfer-matrix. On the lattice and for positive couplings, the 
ground-state energy is given by solutions to the integral equation: 



where ^p{rh) = 1 + Ap/rh, and = 2sin^(p/2), g = b^-g/b^, E = brE, p = bgp and 
m = mb^/b^ One may define a four-particle scattering length a by evaluating the scattering 
amplitude A{p) for four-particle scattering at vanishing external momentum p, and requiring: 

.4-'(0) ^ £ . (12) 

The four-body coupling may then be related to the scattering length by explicit evaluation 
of the inverse scattering amplitude on the lattice. Doing so yields the relation: 

mil 



4vra g g, 



(13) 



c 



where a = a/bg, and gc is obtained by evaluating Eq. 11 at vanishing binding energy. Note 



that for an attractive coupling, by combining Eq. 11 and Eq. 13, one obtains to leading 



order in 1/a (after restoring the lattice spacings) a four particle binding energy: 

This result is very much analogous to that of two particles in three dimensions at large 
positive scattering length, and really follows from dimensional analysis. Particularly, in the 
continuum the only quantity with dimensions of energy that one can construct with available 



scales is Eq. 14 The unitary limit corresponds to tuning the scattering length to infinity. 



or correspondingly, the coupling g to some 0{1) critical value gc- 



^ Throughout this work, I designate dimensionful quantities measured in lattice units with a caret. 
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FIG. 1. Beta function, defined in Eq. 15 plotted as a function of the rescaled bare-coupling 

Note that the lattice mass parameter rh serves as a conversion factor between units of 
length and time, and so one may take hr oc hi providing m is held fixed. The beta function 

= -^7^ , (15) 



can then be computed for the rescaled bare coupling g = g/g^ using Eq. 13 , and by requiring 
that the physical mass and scattering length be invariant under changes of the lattice spacing 
[27] . The result is given by (3g = —g{g — 1) and is plotted in Fig. [ij A continuum theory 
may be defined at the zeros of the beta function; in this case, one finds that there are two 
fixed points: a trivial one at vanishing coupling in the infra-red (IR) corresponding to the 
free theory, and a nontrivial fixed point aX g = gc in the ultra-violet (UV), corresponding to 
the unitary limit. The system is conformal and scale-invariant at both fixed points, and as 
such, no physical scales are available to characterize the system in those limits. 

For this study, I am primarily interested in the nontrivial fixed point located at g = gc- 
Working in the canonical ensemble, every system of fixed total fermion number Q is expected 
to have a zero-energy ground state and a vanishing integrated contact density. This result 
simply follows from the fact that there are no scales in the problem, and therefore all 
dimensionful quantities must vanish. Throughout this work, however, I consider systems of 
fixed fermion number confined to either harmonic trap or a finite box. In such cases, scale 
invariance is explicitly broken by a new length scale that enters into the problem, namely, 
the characteristic size of the system. To unify the discussion for trapped and untrapped 
fermions, I define the characteristic size of the system by Lq = (m/?)"^/^ (trapped) and 
Lq = AL/tc (untrapped). From dimensional analysis considerations, the energy of the system 
must be proportional to the characteristic energy scale u = l/{mLl), and the integrated 
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contact density must be proportional I/Lq. This is true for both the free and interacting 
unitary Fermi gas, although, the constants of proportionality will generally differ in each 
case. Only in the former are the proportionality constants exactly calculable, and in the 
latter case they must be determined nonperturbatively. As discussed earlier, in the unitary 
limit, the energy of the many-body trapped and untrapped systems are given by Eq. [T]and 
Eq. [2} respectively, where the constant of proportionality ^ is undetermined. With the above 
definitions for Lq, the free-gas energies are given by Eq. [4] (trapped) and Eq{Q) = 
(untrapped). At finite volume and in the free theory limit, the integrated contact density 
vanishes both for trapped and untrapped systems. In the former case, this may be seen by 
simply differentiating Eq. |5] with respect to the inverse scattering length and then taking 
a 0". 

In the lattice theory, additional scales appear which also violate the scale-invariance of the 
continuum theory, namely the temporal and spatial lattice spacings. One may quantify the 
lattice discretization errors using the dimensionless parameters = bs/Lo and = brOo = 
rhel. From the latter it is evident that one may independently take the temporal continuum 
limit while holding the spatial lattice spacing fixed by considering the limit m — )■ oo. For 
this study, I fix the anisotropy of the lattice (i.e., fixing rh throughout the study) and then 
extrapolate the characteristic system size (in lattice units) to infinity, or equivalently — )■ 0. 
This procedure allows one to take both spatial and temporal continuum limits (and infinite 
volume limits) simultaneously. Further discussion of the lattice discretization errors and 
their removal may be found in Sec. |Vj 

III. WORLD-LINE REPRESENTATION 

The partition function for the lattice theory is defined in terms of a path-integral over 
fermion fields, weighted by the exponential of the lattice action (often with a chemical poten- 
tial introduced to bias the fermion species numbers toward a desired value). Conventional 
approaches to numerical simulation of the partition function require first reducing the action 
to a fermion bilinear via a Hubbard- Stratonovich transformation (i.e., the introduction of 
non-dynamical auxiliary fields) [281 EH] ■ In the case of the four-body interaction appearing 
in the lattice expression for Eq. |8| this can most easily be achieved using a discrete field 
coupled to ipl^ijjn-er (although there are other equally valid methods as well). One then "in- 
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tegrates out" the fermion degrees of freedom leaving a path-integral over bosonic auxiliary 
degrees of freedom weighted by the exponential of a nonlocal action involving the logarithm 
of a fermion determinant. The resulting effective action is generically complex, rendering 
standard importance sampling techniques which require a probabilistic interpretation for 
the path-integral measure inapplicable. Phase reweighting and other techniques, while in 
principle may be applied to circumvent the problem, are in most cases prohibitively costly 
from the standpoint of computational resources and time due to signal/noise, and other 
problems. 

It was recently demonstrated that a nonrelativistic four- component Fermi gas in one spa- 
tial dimension could be simulated on a lattice free of sign problems by considering alternative 
representations for the partition function [30] . For this study, I use a new path-integral rep- 
resentation, which was inspired by the so-called hopping parameter expansion |3T]. In this 
approach, one may express the partition function for the lattice theory as a path-integral 
over all possible self-avoiding time-directed fermion world-lines. The representation is free 
of sign problems irrespective of population and mass imbalances, making it ideally suited 
for numerical study of few- and many-body four-component fermion systems. Here, I briefly 
summarize the main results of this formulation; for a more in-depth discussion, see [30] . 

For this study, I consider the canonical partition function Z(q) = e~^^^'^'^ for a four- 
component system comprising a fixed fermion number q„ for each species a, and having 
total fermion number given hj Q = 'Yliala- In the world-line representation, the canonical 
partition function is given by the path- integral: 




where C* (g) is the set of all possible self-avoiding loops directed forward in time with a fixed 
winding number q. An example of such a configuration for a single species is provided in 
Fig. |2} For a given configuration c G C*{q), Br{c) represents the total number of time-like 
links associated with the configuration jj and Bs{c) represents the number of space-like links 



^ Note that in general Br{c) = N^q for every c G C*{q). 
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FIG. 2. A possible configuration c G C*(2) for an A^,- = Ns = 6 lattice, reproduced from |30] . 
The set of 11 shaded domains d represent T>{c); seven of those domains have C{d) = 1, three have 
C{d) = 2 and one has C{d) = 3. 



associated with the configuration. T>{c) represents the set of all maximal space-like linear 
domains d of length C{d) formed from the unoccupied sites of c; in Fig. |2| such domains are 
indicated by the shaded gray regions. Finally, the potential term in Eq. 16 is given by 

Va(c) = ^Y1 ^n-|e. ' fi = n + ^e, , (18) 

^ fiSc 

where Va^n = brV^^n- 

In this study, I focus on two physical quantities as a function of the fermion population: 
the energy of the system, E{q) = lim^_>oo F{q), where F{q) is the free-energy of the system, 
and the (integrated) contact density discussed in Sec. [ij Although Monte Carlo simulations 
provide a powerful tool for estimating observable quantities, the free-energy of the system 
itself is generally inaccessible due to the nature of Monte Carlo algorithms. As previously 
noted, however, for fermions in the unitary regime, the only physical length scales in the 
problem are the scattering length a and box size Lq, up to discretization errors. Furthermore, 
the mass parameter is the only quantity available for converting length scales into energy 
scales. One may exploit these observations to gain access to the energy by noting that 
ground-state energy of the system must be proportional to the inverse fermion mass, implying 

m = lim , (19) 
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up to finite lattice discretization errors. 

The energy of tlie system may be determined using tlie Feynman-Hellmann tlieorem, 
whicli may be expressed in tlie path-integral language as: 

dF{q) _ 1 / dS _\ 

n. I ' 
9 



(20) 



(ilogm /3 \(ilogm^ 

where (. . .)g is an expectation value taken with respect to a fixed-charge ensemble associated 
with Note that the full derivative with respect to logm can be written in terms of 

partial derivatives with physical length scales held fixed. Using the relation 

d d f d\ogK\ d f d\ogg\ d , , 

1 a 1 — ™ ) ^11 — . 1 ^11 — ™ ) A 1 — . „ ' ' 



d log m d log m \d log m J d log k, \d log m J ^ d log g ' 
to differentiate the action in Eq. [20] yields three contributions to the energy: 



Eiq)=Tiq) + V{q) + Iiq) , (22) 

in the zero temperature limit. These contributions may be identified as the kinetic (T), 
potential {V) and interaction (/) energies, and correspond to partial differentiation of the 
the action with respect to log(l/m), logn and logg, respectively. In the fermion world-line 
representation, the energy operators are explicitly given by: 

n<i)= ^y/ y ^_iog ""^'^^^y -g.(o\ , (23) 

V{q)=\i^ ' (24) 

tj 

and 



The prefactor appearing in Eq. 25 for the interaction energy operator I{q) at finite scattering 



length may be derived explicitly by differentiating both sides of Eq. 13 with respect to m 



while holding all physical length scales fixed. Doing so yields the useful relation: 

1 / 9 logo \ 1 f dlogQc A 



g \ 9 logm / (jc \ 9 logm 
Taking the temporal continuum limit with bs held fixed (i.e., m oo), one finds that 
Slog^f/^logm — 7- —1, and therefore I{q) — )■ Io{q)- For any finite anisotropy, however, the 
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FIG. 3. Two of the four allowed local constraint-preserving configuration updates (the remaining 
two possibilities are just mirror images of those shown above). 



prefactor remains nontrivial and its inclusion is crucial for obtaining correct continuum limit 
estimates. 

Starting from Eq. [T] and using Eq. 13 , the integrated contact density C (q) (not to be 
confused with the set of closed-loop configurations C* discussed above) may be written as 



C(g) = -{mg) 



,dEiq) 



dg 



(27) 



Taking E{q) as the zero temperature limit of the logarithm of the canonical partition func 
tion, I obtain 

1 



C{q) 
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-{mg) Io{q) 



(28) 



for the contact. One can derive other expressions for the contact by combining Eq. 27 and, 



for example, Eq. [22] as the the expression for the energy rather than the logarithm of the 
partition function. Such expressions are expected to yield the same continuum limit as 



Eq. [28| although estimates based upon such formulae are presumably noisier since they rely 
on the correlations among the various energy observables. For this reason, I only consider 



estimates of the contact based on Eq. 28 



IV. SIMULATION DETAILS 



Monte Carlo simulations were performed for fermions confined to a finite box and a har- 
monic trap. Ensembles were generated using a local updating scheme which preserves the 
constraints placed on the configuration space. Particularly, time-directed links are updated 
one at a time in accordance to the updating rules depicted in Fig. [3j Proposed local updates 
are either accepted or rejected using a Metropolis accept/reject step. Whenever a pro- 
posed update violated the constraints on allowable configurations (such as the self-avoiding 
constraint, or constraints imposed by the lattice boundaries) those proposed updates were 
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TABLE I. Continuum extrapolated observables for up to eight fermions confined to a liarmonic 
trap. 



Q 


qa 


% 




Id 


E/uj (I) 


xVd.o.f. 


E/uj (II) 


xVd.o.f. 


E/io (III) 


xVd.o.f. 


CLo 


xVd.o.f. 


4 


1 


1 


1 


1 


1.008(5) 


1.0 


1.01(1) 


1.1 


1.01(1) 


1.0 


7.13(4) 


0.6 


5 


2 


1 


1 


1 


2.341(6) 


0.1 


2.34(1) 


1.0 


2.34(1) 


0.5 


7.55(5) 


0.7 


6 


3 


1 


1 


1 


4.535(6) 


3.6 


4.50(1) 


1.6 


4.57(1) 


3.0 


8.98(4) 


2.6 




2 


2 


1 


1 


3.610(9) 


0.4 


3.60(2) 


2.8 


3.61(2) 


1.4 


8.25(3) 


2.4 


7 


4 


1 


1 


1 


7.813(6) 


1.1 


7.80(2) 


1.1 


7.83(2) 


1.1 


9.62(4) 


1.2 




3 


2 


1 


1 


5.796(7) 


1.6 


5.79(1) 


1.2 


5.80(2) 


2.0 


9.30(3) 


0.8 




2 


2 


2 


1 


4.719(7) 


0.9 


4.71(1) 


1.5 


4.72(2) 


0.3 


8.87(3) 


1.3 


8 


5 


1 


1 


1 


12.060(6) 


0.9 


12.05(2) 


0.1 


12.07(2) 


0.1 


10.69(2) 


1.2 




4 


2 


1 


1 


9.043(5) 


0.6 


9.03(1) 


0.8 


9.05(1) 


0.6 


10.15(2) 


1.4 




3 


3 


1 


1 


7.964(5) 


0.6 


7.92(1) 


0.6 


8.01(1) 


0.6 


10.34(3) 


1.0 




3 


2 


2 


1 


6.958(8) 


0.9 


6.96(2) 


0.7 


6.96(2) 


1.5 


9.66(3) 


2.0 




2 


2 


2 


2 


4.570(8) 


2.7 


4.55(2) 


6.1 


4.59(2) 


1.3 


19.38(4) 


4.4 



rejected with unit probability. It is known that local updating schemes generically suffer 
from critical slowing, and this updating scheme is by no means any different. However, for 
the lattice volumes and physical parameters explored in this study, the efficiency of the up- 
dating scheme was found to be adequate for achieving percent-level estimates of observables 
with available computational resources. 

All random numbers used in the simulations were generated using Liischer's Ranlux 
pseudo-random number generator with a luxury level equal to one [32]. For this study, 
given the simplicity of the updating scheme, it was found that the random number genera- 
tion was the most time-consuming part of the simulations. The random number generators 
used to generate the configurations for each ensemble were independently seeded so-as to 
yield uncorrelated ensembles. Due to inefficiencies in the Monte Carlo algorithm, however, 
each ensemble involved configurations which were highly correlated in Monte Carlo time. 
Generally, the autocorrelations in each ensemble depend strongly on the simulation param- 
eters considered, and so care was taken to prune the ensembles so-as to eliminate such 
correlations. 
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One interesting feature of the the world hne path-integral representation presented in 



Sec. Ill over conventional approaches that work with fermion determinants, is that in some 
situations there is no computational limitation on the spatial size of the lattice. To see this, 
first note that the computational cost of updating a single configuration by sweeping through 
the lattice scales like PQ. At low temperature, however, the computational cost of a single 
update scales implicitly like the square of whatever length scale in the problem is smallest. 
This is because, roughly speaking, the smallest length scale is what determines the energy 
splittings in the system. So for example, if the only length scale in the problem is the volume 
L, then in order to study the ground state properties of the system, one requires /3 ~ L^. 
If other smaller length scales are present, such as a characteristic trap size Lq of a trapping 
potential or a finite scattering length a, one can then increase the spatial volume arbitrarily 
without increasing the computational cost of the simulation since the energy splittings are 
determined by those other smaller scales. 

The trapped simulations for this study were performed on a finite lattice chosen such that 
L ^ Lq, where Lq corresponds to the trap size. Generally finite volume errors for the trapped 
system depend on the likelihood for the few- or many-body ground state wavefunction to 
lie outside the box [33.j. Given that the ground state wavefunction for trapped unitary 
fermions behaves asymptotically like a harmonic oscillator wave function, one can expect 
finite volume artifacts to be exponentially suppressed in L/Lq. All few- and many-body 
simulations for trapped fermions in this study were performed at spatial lattice volumes 
satisfying L/Lq > 25, and by the scaling arguments above may effectively be regarded 





FIG. 4. Continuum few-body energies (obtained using definition I) and integrated contact densities 
for trapped fermions. The result for the contact determined for q = (2, 2, 2, 2) has been omitted 
from the plot for clarity purposes, but can be found in Table |I} 
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TABLE II. Continuum extrapolated observables for up to 56 fermions confined to a harmonic trap. 
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FIG. 5. Thermodynamic hmit extrapolation of continuum energies and integrated contact densi- 
ties for the trapped Fermi gas. 



as at infinite volume. Although not done so in this study, one may easily monitor the 
configurations as they are updated and verify explicitly whether updates carry fermions 
to the edge of the box when a confining potential is present. The probability for such an 
occurrence during the course of a simulation that has run for a finite amount of time is 
exponentially small in L/Lq. 

Simulations of trapped and untrapped fermions were performed using temporal extents 
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much larger than the expected inverse energy sphttings of the system in order to ensure 
adequate suppression of excited state contamination (or thermal excitations). For trapped 
few- and many-body systems, properties of the Schrodinger algebra imply that the spectrum 
contains a tower of breathing modes, each separated by an amount 2uj [24j. In the case of 
trapped fermions, the temporal extent of the lattice was therefore chosen to satisfy (5uj > 10. 
For untrapped many-body systems, the energy splittings are expected to be of order the 
Fermi energy, given by Ep{Q) = Q'^uj/2, and therefore the temporal extent of the lattice 
was chosen to satisfy l3Ep > 10. Simulations were performed using a single fixed value of 
the lattice mass parameter, m = 1.3, corresponding to a critical coupling ~ 3.7237 

Simulations were performed for multiple values of in order to perform continuum limit 
extrapolations of observables estimated on ensembles of fixed fermion number. Few-body 
trapped ensembles (i.e, Q < 8 for all possible q) were generated for all integer values of 
l/e^ G [3,12], and consisted of approximately 800-50000 uncorrelated configurations after 
thermalizatioi]|^ Note that fewer configurations were generated at smaller as a result of 
increased autocorrelation times due to critical slowing, and also because of the associated 
increase in /3 with Lq. 

Trapped and untrapped many-body ensembles were the same as to those used in [16], 
however, the size of the ensembles have been enlarged, particularly for the untrapped stud- 
ies. Trapped many-body ensembles were generated for Q = 28,32,36,40,44,48,52,56 at 
integer values of l/e^ G [7,20], subject to the constraint Q^^^Cs < 1-0. Ensembles con- 
sisted of approximately 200-1300 uncorrelated configurations with larger ensemble sizes cor- 
responding larger and smaller Q. Untrapped many-body ensembles were generated for 
Q = 32,48,56,64,72,80,88 for equally spaced values of kpbs = Qts G [0.15,7]. Ensembles 
consisted of approximately 600-1600 uncorrelated configurations, again with larger ensemble 
sizes corresponding larger and smaller Q. 

V. ANALYSIS AND RESULTS 

Finite lattice discretization errors may be understood from the view-point of a Symanzik 
action [SU [35] , a continuum description of the lattice theory with lattice spacing dependence 
^ In |16j I erroneously wrote that all simulations consisted of 150-350 configurations; this claim in fact only 
applied to the many-body simulations and not to the few-body cases Q = 4, 5. The latter ensembles were 
essentially the same as those used in this study, and were considerably larger in size. 

17 



carried by the undetermined couplings associated with higher dimension operators. In prin- 
ciple, one should consider the inclusion of all possible local operators consistent with the 
underlying symmetries of the lattice action. Since the continuum theory is scale-invariant, 
dimensional considerations imply that operators of scaling dimension Aq must have associ- 
ated couplings with lattice spacing dependence that scales like b^°~^, where again for a fixed 
physical mass and lattice mass parameter I have used the fact that br ^ b^. Throughout this 
study, I consider dimensionful observables expressed in either units of the trap frequency or 
trap size. Since quantum corrections to the continuum observables typically involve powers 
of coupling constants associated with higher-dimension operators, one may infer the Lq de- 
pendence of such corrections by requiring that the final result be dimensionless. Since Lq 
is the only other dimensionful length scale in the problem besides the lattice spacings, one 
concludes that operators of scaling dimension Aq induce volume dependence scaling like 
gAci-3 ^Qj, dimensionless observables. 

The scaling dimensions of the lowest-dimension few-body operators have been studied 
in detail in [T5] and confirmed numerically in pTBJ. For the theory under investigation one 
concludes that dimensionless observables for a system of fixed total fermion number Q (such 
as the energy measured in units of to or the contact in units of Lq^) must scale as 

OQ(e.) = Oq + 0«e. + O^^^'k'J' + ... , (29) 

where Oq is the physical observable in the continuum limit (independent of m), and Oq^ {j = 
1, 5/3, . . .) are unknown coefficients that depend implicitly on the dimensionless parameter 
rh. The term linear in is the leading spatial lattice spacing error attributed to an untuned 
4 — )• 4 operator with a derivative insertion having scaling dimension Ao/2 = 2. The 
subleading correction is attributed to an parity-odd 5 — ?■ 5 operator with scaling dimension 
Ao/2 = 7/3. Note that the lattice action does not give rise to 2 — )■ 2 or 3 — )■ 3 interactions 
as a result of the point-split nature of the four-body interaction defined in Eq. [TOj 

Continuum limit extrapolations of the dimensionless energies and integrated contact den- 



sities for systems of fixed fermion number were performed using Eq. 29 truncated at order 

5 /3 

es . For trapped systems, three different definitions of the energy were considered: 



E 



T + V + I (/) 

2V {II) 

2(r + j) (///) 
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TABLE III. Continuum extrapolated observables for up to 88 fermions confined to a finite box. 
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The latter two definitions follow from the virial theorem for trapped unitary fermions; the 
first is simply the average of the other two. Although the virial theorem is violated at finite 
e^, in the continuum limit the three definitions should agree. One may either use the three 
definitions to confirm restoration of the virial theorem in the continuum limit for each fixed 
charge system, or one may use the three definitions to gauge the systematic errors in the 
extrapolations. For trapped estimates of the energy, I indicate which energy definition is 
used by the Roman numerals (I), (II) and (III). Energy estimates for untrapped systems use 



definition (I) with y = 0, and all estimates of the contact use Eq. 28 

Fits to observables were performed over a range of values with the maximum varied 
to in order to evaluate the robustness of the fit. Plots of the continuum limit extrapolation 
results for all few- and many-body observables are provided in Appendix [B} Particularly, 
for systems of fixed fermion number I show four plots corresponding to the energy (left) 
and integrated contact density (right). In the upper plots I show the estimated observables 



Oqie^ plotted as a function of es, along with a fit result using Eq. 29 and associated error 
band. The lower plots show continuum extrapolated values Oq obtained for each observable 
as a function of the maximum used in the fit. 

Continuum limit estimates for each trapped few-body observable are tabulated in Table [T] 
along with the corresponding per degree of freedom (d.o.f ) as a measure of the goodness of 
fit. Plots summarizing the fit results are also provided in Fig. |4j For the two cases Q = 4, 5, 
the energies in units of w are known analytically to be one and 7/3 respectively based on the 
operator-state correspondence and knowledge of the scaling dimensions of few-body opera- 

19 



tors \T5]. Continuum extrapolations of the energies yield results statistically consistent with 
the exactly determined values to within 1% and 0.5% statistical errors, respectively. Extrap- 
olation results for the Bertsch parameter and the integrated contact density for trapped and 
untrapped many-body systems are tabulated in Table |TT] and Table III Plots summarizing 
the continuum extrapolation results are shown in Fig. |5] and Fig. |6| 

Thermodynamic limit extrapolations of the continuum observables C,q and Cq (appropri- 
ately normalized) were carried out for the many-body systems. The leading dependence on 
1/Q for these quantities is presently unknown, but is expected to be of the form 



Oq = + 0^p^Q-p + ... . 



(30) 



Following the approach of [16j, I use an ansatz fit function for the Bertsch parameter ex- 
trapolation. For trapped fermions, I fix O^^^ = and determine the parameter O using a 
constant linear least-squares fit over a fit range in which observables appear independent 
of Q (within statistical uncertainties). For untrapped fermions, I use the ansatz p = 1 and 
determine the fit parameters O and 0^p\ Fit results for each case are presented in Table IV 
(trapped) and Table |V] (untrapped) along with the fit range used, and the goodness of fit. 
Fit results and associated error bands are plotted in Fig. |5] and Fig. [6j Within this analysis, 
I obtain the Bertsch parameters ^ = 0.375(1) (I), ^ = 0.374(1) (II) and ^ = 0.375(2) (III) for 
trapped fermions, and C, = 0.370(2) for untrapped fermions. These results are consistent to 
within about two standard deviations, and are also consistent with the analysis of [16]. Sim- 
ilar fits were performed for the contact and yield the values limg^oo CqLq/Q^^'^ = 0.860(1) 
(trapped) and limq^oo CqLq/Q^ = 1.119(5) (untrapped) in the thermodynamic limit. 
Having determined the contacts and integrated contact densities for the trapped and 




0.04 



FIG. 6. Thermodynamic limit extrapolation of continuum energies and integrated contact densi- 
ties for the untrapped Fermi gas. 
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TABLE IV. Thermodynamic limit extrapolated observables for fermions confined to a harmonic 
trap. The three extrapolated values for ^^/^ correspond to measurements made using the observ- 
ables T + V + 1 (upper), 2V (center) and 2{T + V) (lower). Parameters without uncertainties may 
be regarded as fixed during the fitting procedure. 
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TABLE V. Thermodynamic limit extrapolated observables for fermions confined to a finite box. 
Parameters without uncertainties may be regarded as fixed during the fitting procedure. 
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untrapped many-body systems, it is then possible to determine the subleading parameter ( 
appearing in Eq. |5] In the untrapped case the relationship is trivially given by 

CqLq 271 

Plugging in values for the contact obtained from Table [V| yields C, = 0.534(3). A far less 
trivial relation can be derived for the trapped Fermi gas using Thomas-Fermi theory (see 
Appendix |A] for details). Particularly, from Eq. A9 one finds for the trapped case: 



CqLq 8a/2 , 



Plugging in estimates for the trapped contact and Bertsch parameter quoted in Table IV 



yields ( = 0.535(2), which is fully consistent with the untrapped result. Interestingly, one 



may also combine the results of Eq. |3T] and Eq. 32 by equating ( appearing in each formula 
to obtain a third determination of the Bertsch parameter which depends solely the estimates 
of the contact for each system. Doing so yields the value ^ = 0.372(8), which is consistent 
with the other determinations based on estimates of the energy. 
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VI. CONCLUSION 



I have performed lattice Monte Carlo studies of four-component fermion systems confined 
to a finite box and a harmonic trap in one spatial dimension. I presented numerical estimates 
of the energies and integrated contact densities for both few- and many-body systems in the 
unitary limit. The techniques used for this study relied upon a recently developed fermion 
world-line representation for the canonical partition function. The main advantage of this 
representation is that it is free of sign problems for both polarized and unpolarized systems. 
Although not considered here, the representation is also free of sign problems when there is 
a mass-imbalance. 

It was demonstrated in [15] that the unitary four-component gas at zero temperature has 
physical properties that are qualitatively identical to spin- 1/2 fermions at unitarity in three 
dimensions. Numerical studies of the one-dimensional system might therefore provide new 
qualitative and perhaps even quantitative insights into the nature of such nonrelativistic 
conformal field theories. The main findings of this study are: 

1. a less than one percent-level determination of continuum few-body observables for up 
to eight fermions confined to a harmonic trap, providing indirect estimates of the scal- 
ing dimensions of few-body operators based upon the operator state-correspondence; 

2. two independent determinations of the Bertsch parameter ^ to less than one percent 
statistical uncertainties based on estimates of the the continuum ground-state energies 
for trapped and untrapped many-body systems; 

3. a third new determination of the Bertsch parameter to within about two percent sta- 
tistical uncertainties based on estimates of the associated integrated contact densities 
and theoretical input from a calculation based on Thomas-Fermi theory; 

4. two new independent determinations of the parameter ( to within about a half percent 
statistical uncertainties from estimates of the integrated contact densities for trapped 
and untrapped many-body systems; 

5. and finally, verification of the restoration of the virial theorems for all systems of fixed 
fermion number considered. 
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Perhaps the most surprising finding of this study is an apparent numerical equivalence 
of Bertsch parameters for the one- and three-dimensional unitary Fermi gases, the latter 
of which was discussed in the introduction. This observation and its implications were 
originally reported in [16j|. Presently there are no known theoretical arguments for why these 
parameters should be equal. Providing the equality is not by chance, one might naturally 
expect that conformal symmetry and scale invariance plays a crucial role in explaining the 
result. Unfortunately, it unclear whether a similar relationship holds for the contact since 
there is ambiguity in how the scattering length is defined in the one-dimensional system. 

One of the main deficiencies in the analysis presented in Sec. |V] is that although the 
finite volume scaling of dimensionless observables is well-understood from analysis of the 
Symanzik action, presently there is no theoretical understanding of how continuum many- 
body observables depend on the fermion number away from the thermodynamic limit for this 
system. Consequently there is an inherent unquantifiable systematic error associated with 
the thermodynamic limit extrapolations of many-body continuum observables. However, 
the good agreement in ^ and ( obtained from independent untrapped and trapped studies 
provide some confidence that the ansatz fit functions used for the thermodynamic limit 
extrapolation are correct. In the case of the Bertsch parameter, results are further supported 
by a consistent estimate obtained by combining estimates of the contact for each system and 
additional theoretical input. Nevertheless, a theoretical understanding of the scaling with 
1/Q is highly desirable and an obvious place to start for improving the study. 

Secondly, generally speaking, local updating schemes such as the one used in this study 
suffer from critical slowing. It would be interesting to explore whether a worm algorithm, 
or continuous-time Monte Carlo approach could be applied to the nonrelativistic fermion 
world-line formulation in order to improve the efficiency of the simulations. Doing so might 
allow numerical simulations far closer to the continuum and infinite volume limits, and would 
be an important step toward achieving a high-precision (sub-percent level) determination of 
the Bertsch parameter. As previously discussed, such precision determinations could have 
important implications for the three-dimensional unitary Fermi gas as well. 
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Appendix A: Trapped many-body energies from Thomas-Fermi theory 



Here I briefly derive Eq. [2] and its subleading correction in 1/a using a simple density- 
functional theory calculation following [S]. The Thomas- Fermi density functional for har- 
monically trapped fermions is given by 

E^tfIp] = Jdx [S{p{x)) + p{x)v{x)] (Al) 

where 

c 



£{p) = Soip) e 



kpa 



+ ... 



(A2) 



is the energy per unit volume of the system with £o{p) defined in Eq. |3| and v{x) is the 
external harmonic trapping potential defined in Eq. |9} The task is to minimize E!^^p[p] with 
respect to p{x), subject to the constraint that the total number of fermions 

Q = j dxp{x) (A3) 

is held fixed. Introducing a Langrange multiplier p (i.e., a chemical potential) to enforce 



Eq. X3 as a constraint, one finds that the functional is extremized by solutions to 



P 



P = Po 1 - 
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Solving Eq. A4 perturbatively in {pqo) ^, the solution is given by 

.2\ 1/2 



p{x) 



Po 1 



X 



+ 



4C 
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6{xo — \x\) 



(A4) 



(A5) 



(A6) 



up to corrections of order (9(poct)~^- Combining Eq. A3 with Eq. A6, one may relate the 
chemical potential to the charge, finding 

4p 



Q 



+ 



2p 



+ ... 



Inverting this relation yields function of Q, given by: 

1 P2Q^ C 



37raV m ^i/4 



+ 



Plugging Eq. |A6| back into Eq. |A1[ and using Eq. |A5| and Eq. |A8| yields 



^1/4 



(A7) 



(A8) 



(A9) 



for the energy for trapped fermions in the unitary regime, where Eq^^{Q) is given by Eq. [4 
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Appendix B: Continuum extrapolation of observables 



Figures [7 18 show continuum extrapolation results for the energy (upper-left) and inte- 
grated contact density (upper-right) for up to eight fermions confined to a harmonic trap. 
Filled (open) data points in the upper plots indicate data points which have been included 
(excluded) in (from) the fit. Extrapolated results for the energy (lower-left) and integrated 
contact density (lower-right) are shown as a function the maximum eg used for the fit in- 
terval. Boxed data points in the lower plots indicate the extrapolated fit value quoted in 
Table |T] as well as the maximum value use for the fit range in producing the upper plots 
shown in the figure. Dashed lines in the q = (1, 1, 1, 1) and q = (2, 1, 1, 1) plots represent ex- 



act continuum results obtained from [15]. Figures [I9p0| ( 3 1 37) show similar extrapolations 
for up to 56 (88) fermions confined to a harmonic trap (finite box); extrapolated results 
indicated by the boxed data points in the lower figures are provided in Table [ll] (Table 
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FIG. 7. Continuum limit extrapolations for trapped fermions with q = (1, 1, 1, 1). 




FIG. 8. Continuum limit extrapolations for trapped fermions with q = (2, 1, 1, 1). 
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FIG. 9. Continuum limit extrapolations for trapped fermions with q = (3, 1, 1, 1). 
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FIG. 10. Continuum limit extrapolations for trapped fermions with q = (2, 2, 1, 1). 
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FIG. 11. Continuum limit extrapolations for trapped fermions with q = (4, 1, 1, 1). 




FIG. 12. Continuum limit extrapolations for trapped fermions with = (3, 2, 1, 1). 
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FIG. 13. Continuum limit extrapolations for trapped fermions with q = (2, 2, 2, 1). 




FIG. 14. Continuum limit extrapolations for trapped fermions with q = (5, 1, 1, 1). 
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FIG. 15. Continuum limit extrapolations for trapped fermions with q = (4, 2, 1, 1). 
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FIG. 16. Continuum limit extrapolations for trapped fermions with q = (3, 3, 1, 1). 
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q = (3,2.2,1) 
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FIG. 17. Continuum limit extrapolations for trapped fermions with q = (3, 2, 2, 1). 
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FIG. 18. Continuum limit extrapolations for trapped fermions with q = (2, 2, 2, 2). 
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FIG. 19. Continuum limit extrapolations for trapped fermions with q = (3, 3, 3, 3). 




FIG. 20. Continuum limit extrapolations for trapped fermions with q = (4,4,4,4). 
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FIG. 21. Continuum limit extrapolations for trapped fermions with q = (5, 5, 5, 5). 




FIG. 22. Continuum limit extrapolations for trapped fermions with q = (6, 6, 6, 6). 
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FIG. 23. Continuum limit extrapolations for trapped fermions with q = (7, 7, 7, 7). 




FIG. 24. Continuum limit extrapolations for trapped fermions with q = (8, 8, 8, 8). 
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FIG. 25. Continuum limit extrapolations for trapped fermions with q = (9, 9, 9, 9). 




FIG. 26. Continuum limit extrapolations for trapped fermions with q = (10, 10, 10, 10). 
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FIG. 27. Continuum limit extrapolations for trapped fermions with q = (11, 11, 11, 11). 




FIG. 28. Continuum limit extrapolations for trapped fermions with q = (12, 12, 12, 12). 
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FIG. 29. Continuum limit extrapolations for trapped fermions with q = (13, 13, 13, 13). 




FIG. 30. Continuum limit extrapolations for trapped fermions with q = (14, 14, 14, 14). 
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FIG. 31. Continuum limit extrapolations for untrapped fermions with q = (8,8,8,8) 




FIG. 32. Continuum limit extrapolations for untrapped fermions with q = (12, 12, 12, 
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FIG. 33. Continuum limit extrapolations for untrapped fermions with q = (14, 14, 14, 14). 




FIG. 34. Continuum limit extrapolations for untrapped fermions with q = (16, 16, 16, 16). 
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FIG. 35. Continuum limit extrapolations for untrapped fermions with q = (18, 18, 18, 




FIG. 36. Continuum limit extrapolations for untrapped fermions with q = (20, 20, 20, 



42 




O 1.14 



'B.O 0.2 0.4 0.6 0.8 1.0 '''8.0 0.2 0.4 0.6 0.8 1.0 

Q (max.) Q (max.) 



FIG. 37. Continuum limit extrapolations for untrapped fermions with q = (22, 22, 22, 
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